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Abstract 

This  paper  details  the  development  of  a thermoacoustic  mod- 
el and  associated  dynamic  analysis.  The  model  describes  the 
results  obtained  in  a gas  fueled  experimental  combustion  pro- 
gram carried  out  at  UTRC. 

The  contents  of  the  paper  are  (a)  the  development  of  a 
thermoacoustic  model  composed  of  acoustic  and  heat  release 
components,  (b)  the  dynamic  analysis  of  the  resulting  non- 
linear model  using  harmonic  balance  methods  to  compute 
linear  stability  boundaries  and  the  amplitudes  of  oscillations 
and  (c)  the  calibration  of  the  model  to  experimental  data. 

1 Introduction 

This  paper  presents  a thermoacoustic  model  structure  for 
describing  the  dynamics  associated  with  lean  premixed  com- 
bustion instabilities.  The  model  developed  is  used  to  describe 
the  dynamics  of  a single  nozzle  rig  (SNR)  experimental  pro- 


gram carried  out  at  UTRC.  This  paper  also  contains  analysis 
of  this  model  for  both  linear  stability  boundaries  as  well  as 
amplitudes  of  oscillations  subsequent  to  the  loss  of  linear 
stability.  The  calibration  of  model  parameters  using  experi- 
mentally measured  data  is  also  included.  The  purpose  of  the 
modeling  effort  is  to  obtain  a calibrated  model  of  the  single- 
nozzle combustor  that  can  be  used  to  analyze  the  effects  of 
closed  loop  control,  and  most  importantly,  to  ultimately  pre- 
dict the  control  authority  needed  to  reduce  the  amplitude 
of  the  pressure  oscillations  to  an  acceptable  level.  This  pa- 
per represents  a contribution  from  a continuing  program  at 
UTRC  that  develops  combustion  dynamics  and  control  from 
a set  of  scaled  experiments  where  the  scaled  experiments  are 
carefully  constructed  to  capture  the  behavior  of  full  scale  en- 
gine dynamics  and,  importantly,  the  scaled  experiments  may 
be  reflected  back  to  gain  insight  into  the  governing  phenom- 
ena under  investigation. 

The  results  presented  in  the  current  work  both  comple- 
ment as  well  as  extend  related  research  efforts  in  obtain- 
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ing  reduced  order  models  of  combustion  instabilities.  The 
works  of  [5,  7,  8,  11]  are  of  particular  interest.  The  papers 
[5,  7]  concern  the  development  of  coupled  acoustic-heat  re- 
lease models  where  the  heat  release  models  are  focused  on 
unsteady  flame  dynamics.  The  work  of  [11]  is  more  closely 
related  to  the  current  paper  where  the  heat  release  model  is 
focused  on  unsteady  equivalence  ratio  variations.  This  is  also 
true  of  the  work  of  [6]  which  models  afterburner  instabilities 
that  are  driven  by  fluctuating  equivalence  ratio,  however,  the 
coupling  to  the  acoustics  is  different  in  that  work  from  the 
current  paper. 

The  contribution  of  this  paper  is  to  develop  a thermoa- 
coustic model  that  describes  the  single-nozzle  combustor  ex- 
perimental results.  The  methodology  developed  here  consists 
of  describing  the  physics  based  model  components  and  then 
performing  analysis  and  calibration  using  this  framework. 
The  model  framework  consists  of  a feedback  interconnection 
of  an  acoustic  component  driven  by  a heat  release  componen- 
t.  The  nature  of  the  feedback  is  as  follows.  The  fluctuating 
heat  release  component  is  dependent  on  the  fluctuating  e- 
quivalence  ratio  at  a flame  sheet  a distance  away  from  the 
nozzle  exit,  and,  as  it  is  assumed  that  the  fuel  flow  rate  is 
constant  the  fluctuating  equivalence  ratio  is  thus  dependen- 
t on  the  fluctuating  velocity  at  the  nozzle  exit  which  closes 
the  feedback  loop.  The  spatial  distance  of  the  flame  sheet 
from  the  nozzle  exit  can  be  viewed  as  a temporal  time  delay 
(convective  lag)  of  the  nozzle  velocity  to  create  a variation 
in  equivalence  ratio  and  hence  the  structure  of  the  resulting 
model  is  a nonlinear  delay  differential  equation. 

The  resulting  model  of  the  combustion  instability  is  then 
analyzed  using  harmonic  balance  methods  to  compute  linear 
stability  boundaries  and  amplitudes  of  oscillation.  The  mod- 
el contains  an  empirical  gain  factor  and  the  paper  develops 
a methodology  to  calibrate  the  model  to  experimental  data 
using  dynamic  and  stochastic  analysis  tools.  The  metrics  for 
this  calibration  are  (a)  to  match  the  trend  in  amplitude  as 
a function  of  mean  equivalence  ratio,  (b)  to  match  the  fre- 
quencies of  oscillation  and  (c)  to  match  the  power  spectral 
densities  and  probability  function  distributions  of  the  mea- 
sured pressure  time  traces. 

The  organization  of  this  paper  is  as  follows.  Section  2 
presents  the  modeling  framework.  This  includes  both  the 
modeling  of  the  acoustics  as  well  as  the  heat  release  in  the 
combustor.  Section  3 presents  analysis  of  the  model  that 
gives  some  qualitative  insight  into  the  structure  of  the  dy- 
namic behavior  with  special  attention  to  the  self-excited  os- 
cillatory behavior.  Section  4 presents  the  calibration  proce- 
dure that  has  been  used  to  date.  Section  5 summarizes  the 
structure  of  the  models,  the  methodology  that  has  been  de- 
veloped, the  results  from  this  investigation  and  finally  the 
enumeration  of  major  remaining  items  of  investigation. 

2 Thermoacoustic  Model  Structure 

This  section  presents  the  basic  structure  of  the  thermoacous- 
tic model  for  single  nozzle  rig  (SNR).  The  model  is  composed 
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Figure  1:  Single  Nozzle  Rig  Used  in  Experimental  Study 


Figure  2:  Schematic  for  Acoustic  Modeling 


as  a feedback  interconnection  of  acoustic  and  heat  release 
submodels.  The  physical  single  nozzle  rig  is  shown  in  Figure 
1.  The  subsections  below  divide  the  modeling  into  acoustic 
considerations  and  heat  release  considerations. 


2.1  Acoustic  Component 

The  schematic  representation  used  for  modeling  the  single 
nozzle  rig  is  depicted  in  Figure  2.  The  modeling  of  this 
schematic  by  a set  of  coupled  ordinary  differential  equations 
can  be  done  using  momentum  and  mass  conservation  equa- 
tions for  flow  through  the  orifices  and  volumes.  Choking  of 
the  exit  orifice  is  clearly  applicable  as  the  combustor  pressure 
is  roughly  220  psia  and  the  pressure  in  the  exhaust  is  roughly 
20  psia.  The  boundary  condition  at  the  choked  inlet  venturi 
is  taken  to  be  a fixed  mass  flow  rate.  The  equations  for  this 
description  are  as  follows.  Details  may  be  found  in  [15]  for 
the  orifice  relation  and  [4]  for  the  volume  and  choked  orifice 
conditions. 
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The  parameters  are  described  in  Table  1. 

The  resulting  model  is  a nonlinear  differential  equation 
where  the  nonlinearity  occurs  in  the  unsteady  momentum 
equation. 

The  calibration  of  the  model  from  the  mean  values  of  the 
pressure  drops  and  the  air  velocity,  using  mean  flow  rates 
and  cold  flow  (nonreacting)  experimental  data,  is  as  follows. 
The  issue  is  to  choose  the  effective  flow  areas  ( C&A ) of  the 
orifices,  particularly,  the  nozzle  and  the  bypass  leg.  We  can 
use  analytical  formulas  for  steady  states  to  guide  this  choice. 
Steady  state  values  for  pressures  and  velocities  across  the 
system  can  be  computed  using  the  equations  alone: 
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Figure  3:  Combustor  pressure,  nozzle  velocity  and  bypass  leg 
velocity  with  no  combustion 


The  physical  area,  Anozzie,  and  the  effective  flow  area, 
CtfA,  for  the  nozzle  have  been  measured.  The  length  coeffi- 
cients (Lnozzie,  etc)  are  taken  as  physical  parameters  for  the 
physical  length  (but  also  could  be  fit  to  experimentally  mea- 
sured admittance  data  for  the  nozzle  used  in  single  nozzle 
rig  tests).  The  bypass  leg  parameters  are  set  similarly.  The 
discharge  coefficient  (C4)  for  the  bypass  leg  is  taken  as  0.7, 
a typical  value  for  a sharp  edge  orifice.  The  physical  area 
AbyVass  is  adjusted  so  that  21%  of  the  flow  is  through  the 
bypass  leg.  Again,  the  physical  length  of  the  orifice  is  taken 
for  Lbypass-  This  sets  all  parameters  for  the  acoustics.  The 
final  choice  of  parameters  for  the  model  is  given  in  Table  2. 
The  steady  state  values  for  these  parameters  using  (2)  are 
Pc.nmb  “ 228.15,  Ppre  ~ 233.56,  Pplenum  “ 242.07,  Vnozzle  = 
205.69,  Vperf  = 283.86,  Vbypass  = 282.86. 

Calibration  procedure  can  also  include  mean  flow  rates  and 
cold  flow  (nonreacting)  experimental  data.  Specifically  the 
mass  flow  rates  (in  lbm/sec)  predicted  by  the  model  can  be 
compared  to  design  values  for  the  combustor.  The  design 
value  of  mass  flow  rate  is  6.2  lbm/sec  for  the  single  nozzle 
rig  which  compares  with  the  mass  flow  rate  (through  the 
nozzle)  as  5.15  in  the  model  after  calibration.  The  pressure 
drop  across  the  nozzle  was  calculated  as  2.4%. 

A simulation  of  the  equations  of  motion  with  parameter 
values  in  table  2 produces  the  velocity  and  pressure  plots 
shown  in  Figure  3.  The  initial  conditions  were  chosen  to  be 
off  steady  state  to  reveal  a heavily  damped  system, 

2.2  Heat  Release  Component 

This  section  outlines  the  modeling  of  the  heat  release  sub- 
system as  a velocity  sensitive  mechanism. 

The  addition  of  a heat  release  component  to  the  acoustic 
component  considered  above  can  proceed  from  the  equation  for 
mass  conservation  as  in  [10].  The  equation  is 
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where  U is  the  air  velocity  and  p the  density.  The  density 
is  represented  as  a function  of  pressure  P and  entropy  S, 
p = p{P,S).  Then 


dp  _ (dp\  dP  fdp\  35 

dt  \ dP  ) s dt  \ 95 ) p dt 

Then  [10]  the  first  term  can  be  expressed  in  terms  of  the  local 
sound  speed  and  noting  that 


(dp\  =_P_ 

\dSjs  cP 

with  c.p  the  specific  heat  per  unit  mass  at  constant  pressure. 
In  addition  if  there  is  a continuous  volume  distribution  of 
heat  release  at  a rate  of  q per  unit,  mass  then  the  second 
term  is  actually  the  source  term,  that  is, 


The  equation  for  conservation  of  mass  then  becomes 


J_dP 
c2  dt 


-V  • (pU) 


cpT' 


and  the  conversion  to  a volume  average  proceeds  by  inte- 
grating over  the  surface  of  the  control  volume  as  done  in 
[10] . The  main  point  in  the  above  discussion  is  that  the  heat 
release  rate  term  enters  into  the  equation  representing  the 
conservation  of  mass. 

The  assumption  that  is  made  in  this  paper  is  that  the  heat 
release  rate  q is  a function  of  the  velocity  in  the  fuel  nozzle 
q = q(u).  This  modeling  assumption  has  been  discussed 
in  detail  in  [16,  17].  The  basis  of  the  assumption  is  that  the 
unsteady  heat  release  is  a function  of  the  equivalence  ratio  at 
a flame  surface  a distance  l from  the  nozzle  exit,  moreover, 
the  equivalence  ratio  is  only  a function  of  the.  nozzle  exit 
airflow  velocity  if  the  fuel  flow  is  held  fixed.  The  result  is 
that  the  unsteady  heat  release  q(t)  = f (u  (i  — l/i i))  where  u 
is  the  mean  velocity  in  the  nozzle.  The  main  point  to  be  made 
here  is  that  when  a stability  boundary  is  to  be  computed 
for  the  coupled  combustion  system  — the  coupling  of  the 
acoustics  and  the  unsteady  heat  release  — a feedback  loop 
should  exist  between  the  nozzle  velocity  and  the  equation  for 
pressure  fluctuations  in  the  combustor  nozzle. 

The  form  of  the  feedback  can  be  given  more  explicitly. 
The  source  term  F given  in  equation  (1)  that  represents  the 
effects  of  combustion  is  as  follows: 

F(t)  = -*=q  (3) 

Cpl 

with  q being  the  heat  release  rate  per  unit  mass.  The  de- 
velopment in  [17]  can  be  used  to  give  a functional  form  that 
depends  on  the  nozzle  velocity.  A modification  to  the  devel- 
opment in  [17]  is  now  made.  The  modeling  assumption  is 
that  the  heat  release  is  split  into  the  mean  and  fluctuating 
portions  and  that  the  fluctuating  component  depends  only 
on  the  fluctuating  velocity  in  the  nozzle.  Furthermore  the 


heat  release  model  explicitly  splits  the  heat  release  rate  into 
constant  and  fluctuating  part.  This  split  accounts  for  the 
fact  that  in  the  thermoacoustic  model  the  coupling  of  heat 
release  to  acoustics  in  feedback  loop  occurs  through  the  fluc- 
tuating heat  release;  the  corresponding  gain  parameter  for 
the  heat  release  is  then  naturally  introduced  as  multiplying 
the  fluctuating  portion  of  the  heat  release  rate.  The  split 
between  constant  and  fluctuating  heat  release  rate  is  given 
as  follows: 


F{t)  = ~PcAnozzleU&Ha  U (4) 

cpT  \ H(4>)  J 

where  (j>  = is  the  instantaneous  equivalence  ratio 

which  depends  on  fluctuating  nozzle  velocity  ul  = u — u at 
time  t — r,  t is  the  convective  lag. 

In  this  equation  the  heat  release  function  H is  taken  from 
[17]  and  N is  the  gain  parameter  to  be  used  in  comput- 
ing stability  characteristics  of  the  model.  It  is  assumed 
that  the  mean  velocity  u is  equal  to  the  nozzle  velocity  at 
steady  state  conditions,  V^ozzle  (see  equation  (2))  1 . The 
steady  state  conditions  refer  to  the  situation  when  the  rate 
of  change  of  all  variables  is  zero.  Figure  4 presents  typi- 
cal time  traces  for  the  thermoacoustic  model  in  the  insta- 
bility region.  The  time  traces  show  the  limit  cycle  behav- 
ior captured  by  the  model  as  the  system  starts  from  initial 
conditions  and  goes  to  steady  state.  Parameter  values  in 
(4)  are  set  as  cp  = .24  ( BTU / °R  - Ibm),  T = T4  = 3285 
( °R),  pc  = .5807 {lbm/ft3),  MTa  = 1222  {BTU  I Ibm -mix), 
u = 205.69  [ft [sec),  Anozzle  = 0.0431  (/t2). 

3 Thermoacoustic  Model  Analysis: 
Consequences  of  Nonlinearities 

In  this  section  we  analyze  the  dynamics  of  the  thermoacoustic 
model.  We  compute  the  stability  boundary  and  a few  cases 
showing  the  dependence  of  the  amplitude  of  oscillations  on 
parameters  of  the  model.  The  three  varying  parameters  in 
the  model  are  N (gain),  r (convective  delay)  and  <p  (mean 
equivalence  ratio). 

The  harmonic  balance  method  combined  with  parameter 
continuation  (see  Appendix)  is  used  to  compute  the  stability 
diagram  and  the  dependence  of  amplitude  of  the  limit  cycle 
on  model  parameters.  The  N — r stability  diagram  is  shown 
in  Figure  5 and  it  consists  of  several  branches.  Only  one  of 
them  (with  lowest  frequency)  is  actually  computed;  the  oth- 
ers are  obtained  from  it  by  2nk  (k  positive  integer)  shifts  of 

1Varying  $ doesn’t  change  steady  state  values  for  any  of  velocity 
components  but  do  change  all  the  pressures.  Note  that  when  the  system 
is  destabilized  the  unsteady  combustion  can  create  mean  flow  rates 
that  are  different  from  the  steady  ones  because  the  combusting  term  F 
does  not  generally  have  zero  mean.  The  nonlinearity  takes  as  input  an 
unsteady  fluctuating  nozzle  velocity  with  zero  mean  but  will  produce 
as  output  a fluctuating  equivalence  ratio  with  nonzero  mean.  These 
effects  might  be  essential  when  the  instability  is  significant;  in  this 
paper,  however,  the  instability  is  relatively  weak  and  the  drift  of  the  dc 
value  with  gain  can  be  safely  ignored. 
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Figure  4;  Combustion  oscillations  in  thermoacoustic  model 
at  (JV,t,0)  = (.25, 3.23,  .58). 


Figure  6:  Amplitude  diagram  with  varying  gain  (r  = 

3.23,0  = .57). 


StaWftty  diagram  (or  cell  4 SNR  model 


Figure  5:  N — r stability  diagram  for  thermoacoustic  model. 


the  phase  u it  where  w is  the  frequency.  Note  that  as  r in- 
creases u decreases,  which  is  consistent  with  a linear  stability 
analysis  conducted  from  a Nyquist  diagram  of  the  linearized 
system.  The  structure  of  bifurcation  diagram  as  shown  on 
top  panel  Figure  5 has  certain  implications  particularly  there 
will  be  several  instabilities  encountered  as  the  gain  increases. 
This  corresponds  to  multiple  crossings  of  eigenvalues  of  the 
system  into  the  right  half  of  the  complex  plane  which  is  a 
consequence  of  the  infinite  dimensional  nature  of  the  model 
(delay  differential  equation).  In  this  paper  the  focus  is  on  the 
initial  loss  of  linear  stabilty  — the  first  pair  of  complex  conju- 
gate eigenvalues  to  cross  into  the  right  half  plane  — and  the 
ensuing  limit  cycle  behavior  after  this  loss  of  linear  stability. 
The  dynamic  behavior  is  known  as  a Hopf  bifurcation. 

Figures  6,  7 and  8 present  parametric  studies  of  limit  cy- 
cle regime  with  varying  only  one  parameter,  N,  r and  0 
respectively.  The  amplitude  diagram  is  typical  for  super- 
critical Hopf  bifurcation  where  the  amplitude  of  oscillation 
grows  smoothly  (no  jump  behavior)  after  the  loss  of  linear 
stability. 


Figure  7:  Amplitude  diagram  with  varying  delay  ( N = 
.2153,0=  .57). 


AmpMudc  diagram 


Figure  8:  Amplitude  diagram  with  varying  mean  equivalence 
ratio  {N  = .24,  r = 3.23). 
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Figure  9:  Time  traces  for  the  data  used  for  calibration. 


4 Thermoacoustic  Model  Calibra- 
tion to  Experimental  Data 


phi=0.86  phi=O.045  phl=0.624 


Figure  10:  Power  Spectral  Density  for  the  data  used  for  cal- 
ibration. 


This  section  describes  an  approach  for  model  calibration  us- 
ing data  obtained  in  the  experimental  program.  The  data 
used  in  this  section  includes  five  data  points  with  the  mean 
equivalence  ratio  <j>  varying  between  .55  and  .66.  The  mag- 
nitude of  oscillations  reaches  3%  at  the  lowest  4>,  and  the 
frequency  varies  near  220  Hz.  Note  that  the  data  are  quite 
noisy  so  we  should  make  some  assumptions  in  order  to  in- 
terpret the  observed  behavior  as  a limit  cycle  contaminated 
by  noise.  The  main  such  assumption  is  that  oscillations  in 
the  system  are  due  to  the  nonlinear  mechanism  of  a feed- 
back coupling  of  acoustics  and  heat  release  (thermoacoustic 
instability)  which  leads  to  the  emergence  of  a limit  cycle  via 
a Hopf  bifurcation  as  the  equivalence  ratio  becomes  small- 
er. The  equivalence  ratio  decreases  from  the  top  left  to  the 
right  bottom  in  the  figure.  The  most  important  features  of 
Hopf  bifurcation,  for  the  purposes  of  this  study,  are  that  the 
amplitude  grows  as  a square  root  of  the  criticality  parame- 
ter and  the  frequency  varies  linearly  with  the  perturbation 
parameter. 

Time  traces  and  power  spectral  density  of  pressure  signal 
for  chosen  data  points  are  shown  in  Figure  9 and  Figure  10. 
The  identified  values  of  the  magnitude  and  (basic)  frequency 
are  shown  in  Figure  11.  These  qualitatively  agree  with  the 
Hopf  bifurcation  scenario,  but  we  notice  that  the  identifica- 
tion procedure  we  have  used  produces  less  reliable  results, 
especially  in  terms  of  frequency,  for  small  amplitude  oscil- 
lations. This  may  be  explained  by  the  low  signal  to  noise 
ratio.  We  will  consider  two  approaches  to  model  calibration: 
one  for  the  original  thermoacoustic  model  where  no  noise  is 
present  and  another  where  noise  is  included  (into  the  heat 
release  portion)  to  simulate  the  effects  of  turbulence. 


Amplitude  of  oscillations  for  run  48  data  (pts  32-37) 


Figure  11:  Identified  magnitude  and  frequency  of  pressure 
oscillations. 
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Figure  12:  Identified  convective  delay  as  a function  of  equiv- 
alence ratio  (run  48  data,  points  32-34). 

4.1  Model  calibration  without  noise 

The  calibration  procedure  without  noise  proceeds  as  follows. 
We  first  find  a value  of  the  time  delay  as  a function  of  mean 
equivalence  ratio.  The  dependence  of  frequency  on  the  time 
delay  is  used  to  solve  the  inverse  problem:  that  is  knowing 
the  frequency  of  the  limit  cycle  find  the  corresponding  time 
delay.  The  inverse  problem  can  be  easily  solved  graphically 
or  by  interpolating  data  for  the  stability  boundary.  It  is  im- 
portant to  note  that  the  results  of  this  step  do  not  depend 
on  the  particular  choice  of  equivalence  ratio  for  which  the 
stability  boundary  was  computed.  Using  the  fact  that  the 
frequency  of  oscillation  for  available  data,  points  is  known 
the  corresponding  values  of  time  delay  can  be  found.  This 
gives  the  delay  as  a function  of  mean  equivalence  ratio  and 
a parameterization  of  this  function  is  found  by  fitting  a sec- 
ond order  polynomial.  This  allows  the  computation  of  the 
dependence  of  the  amplitude  of  oscillations  on  <j>. 

The  result  of  applying  the  outlined  procedure  to  model  (1) 
and  experimental  data  is  given  in  the  set  of  figures  described 
below.  Table  3 contains  identified  values  of  convective  delay. 
A connection  between  4>  and  r is  described  by  quadratic  de- 
pendence as  illustrated  in  Figure  12.  Using  this  dependence, 
amplitude  curves  have  been  computed  for  several  values  of 
N;  the  results  and  their  comparison  to  the  data  are  shown  in 
Figure  13.  The  comparison  t.o  data  indicates  that  the  data 
fits  best  with  N in  the  range  .25  — .26. 

4.2  Calibration  with  noise  using  statistics  of 
data 

It  is  seen  that  the  pressure  data  has  a large  stochastic  com- 
ponent that  we  will  refer  to  as  the  noise.  Especially  for  the 
low  amplitude  oscillations,  which  occurs  the  higher  mean  e- 
qui valence  ratios,  the  histogram  of  the  time  traces  shows  a 
distinctly  Gaussian  distribution.  Therefore,  we  assume  that 
the  system  is  driven  by  noise.  Note  that  in  [12]  effect  of  ad- 
dition of  noise  to  a combustion  model  was  studied  showing 
good  agreement  between  the  model  and  data.  The  noise  ad- 


Figure 13:  Amplitude  dependence  on  4>  for  the  calibrated 
model  plotted  against  the  experimental  data  (marked  by 
stars).  Three  computed  curves  correspond  (from  left  to  right) 
to  iV  = .27,  .28,  .29. 

dition  to  the  model  was  also  supported  by  results  in  ([2,  3,  9]). 
In  this  work  empirical  models  of  single  nozzle  rig  and  sector 
rig  operating  at  high  equivalence  ratio  conditions  from  forced 
frequency  response  experiment  were  obtained.  A stable  sec- 
ond order  lightly  damped  system  was  identified  from  this  ex- 
periment. The  model  was  augmented  by  addition  of  a white 
noise  input  to  match  the  open-loop  PSD  of  the  combustor 
pressure  signal.  Excellent  matches  of  pressure  PSD  from  the 
model  simulation  and  from  experiment  for  controlled  and  un- 
controlled system  were  obtained.  A general  theory  for  com- 
parison of  models  of  systems  with  complex  non-equilibrium 
attractors  and  stochastically  driven  systems  with  applica- 
tions to  model  validation  and  parameter  identification  is  laid 
out  in  [14].  This  paper  also  contains  an  example  of  parameter 
identification  (delay  and  noise  standard  deviation)  of  a semi- 
empirical  combustion  model  of  the  controlled  single-nozzle 
rig  operating  at  high  mean  equivalence  ratio  using  statistics 
of  data. 

In  order  to  deal  with  this  stochastic  input  explicitly,  a 
random  fluctuation  is  added  to  the  velocity  in  the  nozzle. 
This  random  component  models  turbulence. 

If  no  noise  were  present,  as  mentioned  earlier  in  the  paper, 
the  equivalence  ratio  would  be  <p(t)  — <j>u/u(t  — T).  The  noise 
is  introduced  into  it  in  the  following  manner 

<t>(t)  = 4>ul  ( u(t  - r)  + Nnur(t))  (5) 

where  r is  a normally  distributed  random  number  with 
variance  1.  The  parameter  Nn  determines  the  level  of  input 
noise  entering  the  system  (as  a percentage  of  the  mean  nozzle 
velocity). 

By  comparing  statistics  of  the  data  with  the  simulated  out- 
put of  the  model,  a fit  can  be  done  on  the  model  parameters. 
With  noise  addition,  the  parameters  of  the  model  are  gain 
N,  time  delay  r and  noise  level  N„.  Statistics  compared  are 
the  Power  Spectral  Density  (PSD)  and  the  Probability  Den- 
sity Function  (PDF).  Since  the  model  can  only  describe  the 
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Figure  14:  The  combined  heat  release  function. 


bulk  acoustic  mode,  and  the  data  shows  strong  component 
around  750Hz  corresponding  to  a longitudinal  mode,  the  da- 
ta is  filtered  around  the  bulk  mode  frequency  before  com- 
parison. Six  data  sets  are  available  for  six  different  mean 
equivalence  ratios,  and  for  each  of  these  data  points,  values 
of  these  parameters  are  determined  which  make  the  model 
reproduce  the  characteristics  of  the  data  closely  (PSD  and 
PDF).  In  particular,  six  values  of  N and  r,  and  one  N„  are 
determined.  These  are  shown  in  Table  4. 

Since  the  heat  release  rate  given  by  F(t)  in  equation  4 has 
N and  <j>  as  parameters,  this  would  mean  a different  func- 
tion is  required  to  fit  the  data  for  every  data  point.  For 
this  reason,  a single  heat  release  function  was  constructed 
with  the  help  of  these  six  identified  functions  so  that  when 
the  model  was  simulated  with  it,  the  model  output  again 
matched  the  data  statistics  as  well  as  it  did  with  the  individ- 
ual iV-dependent  functions.  Figure  14  shows  the  heat  release 
function.  Figure  15  shows  the  comparison  of  PSD  and  PDF 
estimates  between  the  data  and  the  model  output  with  this 
heat  release  function  ( Nn  and  r values  used  were  the  ones 
given  in  Table  4).  These  plots  show  that  the  model  with 
the  identified  heat  release  function  can  reproduce  the  chosen 
statistics  of  the  data  reasonably  well.  Figure  16  additionally 
confirms  this  by  showing  the  comparison  of  rms  as  a function 
of  i from  data  and  the  model. 


5 Conclusions 

This  paper  has  presented  a thermoacoustic  modeling  frame- 
work that  describes  the  results  obtained  in  an  experimental 
program  conducted  at  UTRC. 

The  paper  details  the  overall  modeling  approach  including 
the  development  of  physics  based  acoustics  and  heat  release 
components.  The  feedback  interconnection  of  acoustics  and 
heat  release  — the  thermoacoustic  model  — has  been  an- 
alyzed using  a harmonic  balance  method.  This  analysis  is 
used  to  calibrate  the  model  to  experimental  data,  especially, 
to  match  amplitudes  and  frequencies  of  oscillations  seen  in 


Figure  15:  Comparison  of  data  characteristics  with  that  of 
model  output  with  best-fit  heat  release  function  (10%  noise). 
Left  column:  PSD’s;  right  column:  PDF’s. 


Pressure  amplitude  (rms)  (or  data  and  best-fit  model 


Figure  16:  Pressure  amplitude  (rms)  for  run  48  data  and 
best-fit  model  with  noise  (noise  level  is  the  same  at  al- 
1 points.) 
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experiments.  The  matching  is  accomplished  by  selection  of 
the  empirical  portions  of  the  model  — the  time  delay  and  the 
heat  release  gain.  The  confirmation  or  match  of  the  model  to 
trends  in  the  amplitude  of  oscillation  as  a function  of  mean 
equivalence  ratio  forms  the  validation  of  the  model. 

This  paper  has  developed  a methodology  for  both  mod- 
eling and  analysis  that  captures  experimental  data.  This 
methodology  consists  of  the  following  steps: 

1.  The  development  of  an  acoustic  model  based  on  con- 
servation laws.  The  acoustic  model  captures  the  bulk 
motion  in  the  combustion  chamber.  The  model  uses  re- 
alistic boundary  conditions  for  the  premixing  nozzle  and 
the  exit  orifice.  The  calibration  of  the  acoustic  model  is 
done  using  steady  state  rig  measurements  of  pressures 
and  mass  flow  rates  and  the  model  predicts  the  dynam- 
ic behavior  including  acoustic  resonance  frequencies  and 
damping. 

2.  The  development  of  a heat  release  model  based  main- 
ly on  first  principles.  The  heat  release  model  captures 
the  mean  heat  release  as  well  as  the  fluctuating  com- 
ponent. The  fluctuating  component  is  a function  of  a 
convective  time  delay  and  an  empirical  gain  that  both 
need  to  be  identified  from  experimental  data.  The  mean 
heat  release  is  a function  of  the  fuel  mass  flow  rate  and 
equivalence  ratio  and  can  be  calculated  separately. 

3.  The  calibration  of  the  empirical  portions  of  the  model  to 
match  experimental  data.  The  paper  presents  two  meth- 
ods for  calibration  of  the  model.  In  the  first  method  the 
time  delay  is  estimated  from  experimental  data  to  match 
the  frequency  of  oscillation  and  is  expressed  as  a func- 
tion of  mean  equivalence  ratio.  The  heat  release  gain  is 
chosen  so  that  the  amplitudes  of  oscillation  of  the  model 
matches  the  experimental  data  at  one  value  of  equiva- 
lence ratio.  The  model  is  then  used  to  predict  the  trend 
of  amplitudes  with  respect  to  equivalence  ratio  and  this 
validates  the  model  against  experimental  data.  The  val- 
idation shows  good  agreement.  In  the  second  method 
a stochastic  fluctuation  in  the  fuel  nozzle  velocity  is  in- 
troduced and  the  PSD  and  PDF  properties  of  the  data 
are  matched  through  adjustment  of  the  noise  and  gain 
values.  The  validation  shows  better  agreement  with  the 
data  than  the  first  method  particularly  in  the  low  am- 
plitude portion  of  the  pressure  data  where  the  system  is 
expressed  as  a noise  driven  stable  linear  system. 

The  paper  enables  a number  of  subsequent  investigations. 
Some  of  the  immediate  ones  are  as  follows: 

1 . Use  of  closed  loop  data  for  model  calibration  and  vali- 
dation. 

2.  Prediction  of  control  authority  to  match  experimental 
data  in  a quantitative  manner. 

3.  Use  of  forced  response  data  to  validate  the  model  against 
experimental  data. 


4.  The  calibrated  model  allows  the  use  of  control  theory 
to  understand  the  fundamental  limits  of  controlled  per- 
formance. This  is  a standard  use  of  control  theory  as 
developed  in  [18]  and  applied  to  combustion  processes 
in  [2,  3],  This  study  of  the  limits  of  controlled  behav- 
ior allows  the  central  parameters  to  be  understood  and 
controllers  to  be  evaluated  against  absolute  metrics  and 
not  relative  advantages  of  control  design  methods. 

All  of  these  investigations  are  expected  to  increase  the  fi- 
delity of  the  model  by  examining  possible  deficiencies  and 
through  the  incorporation  of  alternate  parameter  selection  or 
physical  mechanisms  more  accurately  describe  a wider  range 
of  experimental  conditions  than  the  open  loop  matching  that 
is  contained  in  this  paper. 
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Appendix:  Harmonic  balance 

method 

The  purpose  of  this  appendix  is  to  provide  a computational 
framework,  using  harmonic  balance  method,  for  tracing  sta- 
bility boundaries  and  tracking  periodic  orbits  as  parameters 
in  the  model  vary.  We  assume  the  model  to  be  in  the  form 
of  differential-delay  equation: 

x = F (x(t),x(t  — r),  a)  (6) 

Here  x € Rn  denotes  the  state  vector,  r is  the  system  delay, 
and  a is  a scalar  parameter. 

Let  x — X(t)  be  a periodic  solution  of  the  system  (6), 
X(t  + T)  = X{t),  where  T is  an  unknown  period.  Assuming 
that  x(t)  can  be  approximated  by  the  truncated  Fourier  series 
with  M harmonics  then 

M 

Xm  {t)  = xo  + ( Pk  sm(2nkuit)  + cos(27rku;t))  (7) 

*=1 

where  w = 1/T,  x0,pk,qk  G Rn-  Substituting  (7)  into  (6) 
and  requiring  the  truncated  Fourier  series  of  both  sides  of 
(6)  be  equal  yields  the  following  set  of  M-th  order  harmonic 
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balance  equations: 


where 


F0 


FSh 


FCk 


F0 

= o, 

FSk  - 2n kuiQk 

= 0,(& 

FCk  + 2nbupk 

= 0 ,(k 

(pi  1 9i ) 

= 0. 

= 1,  ...M), 

= 1 

(8) 


u f F(XM(t),XM(t  ~ r),a)dt, 

Jo 

rl/u 

2 ui  / F(Xm (t),  Xm (t  — r),  ci)  sin(27r/c£jt)dt, 

Jo 

/•1/u 

2u  j F(Xtt(t),  X\i(t  — r),  a)  cos(2Trkut)dt, 

Jo 


(9) 


Note  that  the  last  equation  in  (8)  fixes  the  phase  thus 
removing  the  singularity  due  to  translation  invariance  of 
the  periodic  solution.  The  algebraic  system  (8)  contains 
(2 M 4-  l)n+  1 equations  and  the  same  number  of  unknowns 
(vectors  xq ,Pk,Qk,  (k  = 1,  ...M)  and  frequency  u>).  Taking  r 
as  another  unknown  yields  a standard  continuation  problem 

[1]  since  the  solution  set  of  (8)  will  generally  be  a curve.  It 
determines  the  evolution  of  a periodic  solution  as  r varies. 
Note  that  instead  of  r we  may  choose  a as  an  active  param- 
eter and  analyze  the  dependence  of  the  periodic  solution, 
particularly,  its  amplitude  and  frequency,  on  a. 

For  the  purposes  of  linear  stability  analysis  linearize  equa- 
tion (6)  at  the  equilibrium  point  (which  is  to  be  determined) 
and  find  critical  value  of  delay  r and  a such  that  linearized 
equation  has  periodic  solutions.  Since  these  solutions  are 
pure  harmonics,  we  use  (8)  with  M — 1.  We  also  have  to  fix 
the  amplitude  of  oscillations  (which  because  the  system  is  a 
linear  system  can  be  arbitrary,  moreover  is  set  by  the  initial 
conditions  of  the  defining  system  of  equations)  which  we  do 
by  adding  the  constraint  (f>i,Pi)(9i,9i)  — (pi,9i)2  = 1-  In 
more  detail,  the  defining  system  for  stability  calculations  is 
given  as  follows. 

The  linearized  system  (6)  at  the  equilibrium  xo 
(F(xo,  Xq.  a)  = 0)  has  the  form 

^7^  = Z(xQ,a)u(t)  4-  W(xQ,a)u(t  — r)  (10) 

at 

where  Z and  W are  tixn  matrices.  The  defining  system  now 
reads: 


F(x  0)  = 0, 

Zp i 4-  NW ( pi  cos(27twt)  4-  qi  sin(2jrwr))  + 27r<q  = 0, 

Zqi  4-  NW(—pi  sin(2jT6L)r)  4-  qi  cos(2irwr))  — 2irpi  = 0, 

(Pi  j 9i ) = 0, 

(pi,Pi)(9i,9i)  = 1(H) 


This  system  contains  3n4-2  equations  and  the  total  number  of 
unknowns  (xo , pi , 9i , w)  is  3n  + 1 . Adding  r and  a as  another 
unknowns  yields  a continuation  problem  for  computing  the 
onset  of  (linear)  instability.  The  projection  of  the  solution 
curve  onto  (a,  r)  plane  defines  the  a-r  stability  diagram.  The 


evolution  of  frequency  along  the  stability  boundary  can  be 
seen  from  projections  of  the  same  curve  onto  (r,  c j)  or  (a,  u) 
plane. 

We  now  turn  to  the  problem  of  computing  limit  cycles  of 
small  amplitude  near  the  stability  boundary.  Limit  cycles 
near  Hopf  bifurcation  [13]  can  be  approximated  with  two 
harmonics  (M  — 2).  The  corresponding  algebraic  system  of 
harmonic  balance  equations  contains  5ra  + 1 equations  and 
unknowns.  Adding  r (or  a)  to  the  set  of  unknowns  yield- 
s a continuation  problem.  The  solution  of  these  equations 
determines  the  dependence  of  the  amplitude  and  frequency 
of  oscillations  on  r.  As  we  move  away  from  the  onset  of 
instability  into  the  unstable  region,  adding  more  harmonics 
might  be  necessary  to  represent  the  limit  cycle  (particularly 
when  saturation  effects  of  nonlinearity  become  important) 
with  sufficient  accuracy. 
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Parameter 
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density  at  the  nozzle  exit, 
(lbrn/ft3) 

Lnozzlei  Lperf,  Lbypa8S 

effective  nozzle  (perforated 
plate,  bypass  leg)  length, 
(ft) 

^nozzle  i kper/i  & bypass 

\/C'd  for  nozzle  (perforated 
plate,  bypass  leg)  flow,  (di- 
mensionless) 

Vnozzle  ? ^ perf  5 V bypass 

velocity  at  nozzle  (perforat- 
ed plate,  bypass  leg)  exit, 
(ft/sec) 

P combi  ^plenum  i Ppre 

pressure  in  the  combustor 
(plenum,  prediffuser),  (psia) 

C combi  C 'plenum  » C pre 

capacitance  of  combustor 
(plenum,  prediffuser)  vol- 
ume, (ft  — sec?) 

Ae 

physical  area  of  exit  cross 
section,  (ft2) 

K 

constant  in  choked  flow 
equation,  (V  °R/sec) 

cd 

discharge  coefficient,  (di- 
mensionless) 

T4 

exit  temperature  in  combus- 
tor, ( °R) 

rhin 

upstream  mass  flow,  (lb- 
m/sec) 

F 

mass  flow  addition  due  to 
heat  release(lbm/sec) 

Table  1:  Description  of  Model  Parameters 


Parameter 

Value 

Comments 

P3 

.5807 

lbm/ft3 

L nozzle 

.65 

ft 

L bypass 

.95 

ft 

Lptrf 

.042 

ft 

knozzle 

2.0408 

dimensionless 

R( bypass 

2.7778 

dimensionless 

kperf 

1.6866 

dimensionless 

■A-nozzle 

.0431 

ft2 

-A-bypa.88 

.0083 

ft2 

Aper  f 

.0313 

ft2 

Ccomb 

2.7455  x 10~4 

ft  — sec2 

Cpre 

1.8122  x 10“* 

ft  — sec2 

^ plenum 

4.1853  x 10“s 

ft  — sec2 

6.52 

lbm/sec 

K 

.532 

R/sec 

cd 

.98 

dimensionless 

Ae 

.0218 

ft2 

t4 

2825+460 

° R 

Table  2:  Parameter  Values  for  Acoustic  Portion  of  Model 
(!)• 
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Mean 

equivalence 

ratio 

Amplitude 
(percent  to 
the  mean) 

Frequency 

(Hz) 

Identified 
delay  (ms) 

0.5520 

3.25 

216 

3.23 

0.5750 

2.3 

221 

3.11 

0.5990 

1.149 

225 

3.02 

0.6240 

0.5075 

222 

N/A 

0.6450 

0.3410 

223.7 

N/A 

0.6600 

0.2795 

223.5 

N/A 

Mean 

equivalence 

ratio 

Best-fit  N 

Best-fit  r 
(ms) 

Nn  (noise 
level) 

0.66 

0.278 

3.0 

0.1 

0.645 

0.272 

3.0 

0.1 

0.624 

0.259 

3.0 

0.1 

0.599 

0.244 

3.0 

0.1 

0.575 

0.22 

3.08 

0.1 

0.552 

0.19 

3.20 

0.1 

Table  3:  Data  from  run  48  and  identified  value  of  delay  Table  4:  Best  fit  values  of  N , r and  Nn  to  data  from  run  48. 


PAPER  -30,  C.  Jacobson 

Question  (M.  Huth,  Germany) 

How  can  models  developed  for  a single  burner  test  rig  be  used  in  an  annular  combustion 
arrangement?  How  are  the  circumferential  modes  described? 

Reply 

The  UTRC  approach  has  been  to  develop  combustion  dynamics  models  by  reducing  the 
number  of  fuel  injectors,  but  keeping  the  environment  (pressure,  temperature)  realistic 
relative  to  the  engine  operating  conditions.  Hence  the  heat  release  models  are  realistic 
and  the  acoustic  model  may  be  changed  to  model  different  situations  (longitudinal,  tan- 
gential modes). 

Question  (S.  Candel,  France) 

Can  you  explain  the  meaning  of  the  gain  parameter  and  is  this  gain  parameter  a function 
of  frequency?  Also,  in  many  cases  the  response  of  a confined  flame  is  strongest  in  re- 
gions where  the  gain  of  the  transfer  function  is  not  very  large  but  where  the  phase  is 
large. 

Reply 

The  dynamic  model  developed  uses  physics-based  models.  However,  these  models  were 
very  simplified.  In  the  procedure  described  in  the  paper,  a free  parameter  (“gain”)  is  in- 
terested to  capture  some  neglected  heat  release  characteristics.  This  gain  parameter  is 
expected  to  be  0(1)  if  the  model  is  accurate  (which  it  was  found  to  be  when  N « 0.2  in 
the  calibration  procedure). 

The  heat  release  model  was  derived  based  on  the  work  of  Proscia  and  Peracchio.  In  this 
model,  both  equivalence  ratio  and  flame  sheet  dynamics  are  considered  and,  in  this  full 
model,  there  will  be  a frequency  dependence  (because  of  the  flame  sheet  dynamics).  This 
dependence  was  not  found  to  be  significant  for  the  behavior  being  captured  by  the  model 
and  hence  was  deleted  to  have  the  simplest  model  available. 


